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ABSTRACT 

Young star clusters like R136 in the Large Magellanic Cloud and NGC 3603, Westerlund 1, and 
2 in the Milky Way are dynamically more evolved than expected based on their current relaxation 
times. In particular, the combination of a high degree of mass segregation, a relatively low central 
density, and the large number of massive runaway stars in their vicinity are hard to explain with the 
monolithic formation of these clusters. Young star clusters can achieve such a mature dynamical state 
if they formed through the mergers of a number of less massive clusters. The shorter relaxation times 
of less massive clusters cause them to dynamically evolve further by the time they merge, and the 
merger product preserves the memory of the dynamical evolution of its constituent clusters. With a 
series of A^-body simulations, we study the dynamical evolution of single massive clusters and those 
that are assembled through merging smaller clusters together. We find that the formation of massive 
star clusters through the mergers of smaller clusters can reproduce the currently observed spatial 
distribution of massive stars, the density, and the characteristics (number and mass distribution) of 
the stars ejected as runaways from young dense clusters. We therefore conclude that these clusters 
and possibly other young massive star clusters formed through the mergers of smaller clusters. 
Subject headings: methods: numerical - open clusters and associations: individual (R136, NGC 3603, 
Westerlund 1, Westerlund 2) - galaxies: star clusters 



1. INTRODUCTION 

Young ( < 4 Myr) and dense {p > lO'^Mo/pc^) 
star clusters, NGC 3603, Westerlund 1 and 2 in the 
Milky Way Galaxy and R136 in the Large Magellanic 
Cloud (LMC), pose a number of interesting evolution- 
ary problems on the formation and early evolution of 
star clusters. Their relatively small number of stars and 
their young ages make them ideal modeling targets (see 
iPortegies Zwart et al.|[2QTQl ). because they can be simu- 
lated on a star-by- star basis for their entire lifetime. In 
addition, the observational data on these clusters have 
excellent quality, which makes them ideally suited for 
comparison with the results of the numerical simulations. 

From an observational point of view, each of 
these clusters shows a remarkably mature dynamical 
state compared to their ages (Harayama et al. 2008; 
[Andersen et al.l l2009l : lAscenso et al.ll2007l : IGennaro et al.l 
Hpil); the degree of mass segregation suggests that either 
these clusters are born with a certain degree of mass seg- 
regation or that it develops more quickly than expected 
by standard relaxation theories. 

In particular, the LMC cluster, R136, shows 
characteristics of a post-core-collapse star cluster 
(jMackev fc Gilmord[2003l : iFuiii fc Portegies Zwartll201lL 
hereafter FPZ2011) . The cuspy density profile of R136 
(jMackev fc GilmQrdl2003f ) and the large number of mas- 
sive runaway stars which seem to have escape from R136, 
such as VE TS 682 (Bestenlehner et al. 201 1 |, 30 Dor 016 
(lEvans et al .1120101) and severa l OB stars (iBrandl et al.l 
120071 : iGvaramadze "ind1[20T0l ). are effective signatures 
of core collapse. The mass and number distribution of 



the runaway stars ejected from R136 can be explained 
by three-body interact ions between a hard binary and a 
single star (FPZ201 1; IGvaramadze fc GualandrisI 120111 : 
iBaneriee et ani2012f ). Such a "bully" binary (BB) that 
formed in the gravothermal collapse of the cluster kicks 
out its surrounding stars and produces runaway stars. In 
FPZ2011, we demonstrated that the number of runaway 
stars that one BB can produce does not depend on the 
global parameters of the cluster, and accounts for 21 
in total. The fraction of runaway stars to the whole clus- 
ter, /run, is therefore larger in clusters with fewer stars. 
The observed fraction of runaway stars is satisfactorily 
explained by the cluster having experienced core collapse 
within 1 Myr since its birth (see FPZ2011). 

Such an early core collapse can be initiated suffi- 
ciently quickly if the initial density of the cluster core 
pc ^ lO^M0pc~^, which is considerably higher than the 
currently density observed in R136: pc ~ 5 x IO^MqPc"^ 
(|Mackev fc Gilmord 120031) . With a current half-mass 
relaxation time of trh ^ 100 Myr (jMackev fc Gilmord 
l2003f ). it is unlikely that R136 evolved to a state of core 
collapse within 1 Myr and subsequently evolved towards 
the less dense state it is in today. 

One way to speed up the initial dynamical evolu- 
tion towards a state of core collapse is by forming 
the cluster from the hierarchical mer ging of several 
smaller clusters (jMcMihan et al.l I2007D . The shorter 
relaxation timescales of these smaller clusters cause 
them to evolve dynamically more quickly than one sin- 
gle large cl uster (lA arseth fc Hills 1972; McMill an et al.l 
120071 : Moeckel fc Bonneh 2009; Yu et al.l 120111 ^ A 
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merger between these dynamically evolved clusters pre- 
serves the memory of their past dynamical evolution and 
impri nts this onto the more massive post-merger clusters 
(McM illan et"all l2007l ) . 
Sub-clustered star formation seems to be common 
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for young stars and protostars embedded in molecu- 
lar clouds (iGutermuth et all I2QQ9I: iLada Lad"al I2QQ3I : 
iPortegies Zwart et alJ l2QlQf ). Numerical simulations of 
star formation in turbulent molecular clouds also sup- 
port star cluster formati on via hierarchical mergers 
(iBoimen et al. 2003 2QQ4 [20TTI : [Girichi dis et al.. 2D1^ 
iKruijssen et al.fl2012f ). Not only for young clusters in 
the Galactic disk, the merger scenario is suggested for 
extended old globular clusters in the Milky Way - like 
NGC2419 (jBriins &: KrouDall2011l :lBruns et al.| [20TTf ) and 
massive star clusters w hich are born during galaxy merg- 
ers ([Saitoh et al.ll20H"h . Thus, the formation of star clus- 
ters via hierarchical merging of sub-clusters might be 
very common for any kinds of star clusters. 

The hierarchical merging of small sub-star clusters 
solves two problems in our understanding of R136 and 
other clusters with a similar age, density, mass, and mass 
function: (1) It allows for the cluster with a long re- 
laxation time to experience core collapse well within the 
time frame for these sub-clusters, and (2) because each of 
the sub-clusters experiences an individual core collapse 
before they merge, the number of runaway stars pro- 
duced by dynamical sling shots is considerably higher. 
Since the number of runaways is independent of clus- 
ter mass (FPZ2011), the expected number of 0-type 
runaway stars is only one or two for R136 if it formed 
as a single cluster, but at least three were observed 
(|Gvaramadze et al.ll20Tol : lEvans et"aIH 2010). It is there- 
fore surprising that R136 has produced more runaways 
than anticipated by this theory if it was born as a sin- 
gle cluster. The hierarchical merging formation of R136 
solves this conundrum. 

The hierarchical merging of star clusters is also prefer- 
able for the formation of very massive (> 150 Mc^ ) stars 
found in R136 and NGC 3603 (Crow ther et ldTl2010). 
Runaway collisions of stars during core collapse eas- 
i ly form such very niassive stars in the clu ster core 
(|Portegies Zwart et al.l 119991 : iFuiii et al.l I2009D . These 
massive stars stay in the cluster center due to mass seg- 
regation and form BBs. In this way, each sub-cluster 
that experiences core collapse will lead to the forma- 
tion of one BB. After the sub-clusters merged, the mas- 
sive BBs interact with each other in the core of the 
merged cluster. The binary-binary i nteractions can form . 
a massive escaping; binary like R1 45 (jSchnurr et al .1120091 : 
iGvaramadze fc Gualan dris 2011). 

We also argue that Westerlund 1 and 2 are still 
in the hierarchical merging process, which is demon- 
strated by their clumpy appearance. These clusters 
have not completed their merging processes, which 
has profound consequences for their current degree of 
mass segregation a nd the distri but ion of their stars 
(iMengel &: Taccon i-GarmanI [20091 : iGennaro et al.l 120111 : 
iRoman-Lopes et a l. 2011). 

In this paper, we demonstrate that merging star clus- 
ters can explain the dynamically evolved characteristics 
of young dense clusters in the Milky Way and the LMC. 
The observed distributions of massive stars including 
stars with an initial mass of ^ ISOM©, such as those 
found in R136 and NGC 3603, which are consistent with 
our results. Furthermore, the runaway stars around R136 
and Westerlund 2 are consistent with our merger models. 
The core collapse in each small sub-cluster before merger 
forms a massive hard binary in its core. These massive 



binaries in the sub-clusters produce a sufficient amount 
of runaway stars, and then they are ionized or ejected 
from the cluster after their hosts merged. We demon- 
strate that such massive stars exhibit a spatial distribu- 
tion similar to those observed in the region of R136 and 
NGC 3603. These dynamical processes result in a rela- 
tively low density of the cluster compared to their initial 
density or the density during the core collapse; the core 
density of these merger remnants in our simulations is 
consistent with those of observed young clusters. 

2. METHODS 

We performed a series of A/'-body simulations of sin- 
gle and multiple star clusters that merge to form a large 
single cluster. For merger cases, we adopted a mass of 
6.3 X lO^Mo and a King model (King 1966) with the 
dimensionless central potential Wo = 2 for each sub- 
cluster. The initial mass function (IMF ) of stars were 
drawn from the Salpeter mass function ()Salpeterlll955f ) 
between 1 and 100 Mq. We imposed that each clus- 
ter contains at least one star with m > SOMq. The 
initial half-mass radius and core density of our simu- 
lated clusters are rh — 0.1 pc and pc — 2 x 10^ Mq 
pc~^, respectively. This core density is similar to the 
central density of sub-clusters formed in star formation 
simulations in a turbulent molecular cloud ([Bonnell et al.l 
12004 ). With these initial conditions our model clusters 
have = 2048 (2k) stars and a half-mass relaxation 
time of trh — 0.37 Myr. We call this model single-2k. 
The parameters are summarized in table [H 

We initially distribute 4 or 8 of these sub-clusters 
(single-2k) randomly in a volume with a radius of Tmax 
and with zero velocity. We varied rmax from 1 pc to 6pc. 
The clusters merge within a few Myr to a single cluster. 
The total masses of the merger remnant correspond to 
the mass of NGC 3603 for 4-cluster mergers and R136 
for the 8-cluster merger cases. We summarize the initial 
conditions in table [2l 

For comparison, we also performed simulations of sin- 
gle clusters, whose masses are identical to R136. For 
this single cluster model, we adopted a King model for 
the initial density profile with Wo = 6; we call this model 
single- w6. The core density of model single- w6 is com- 
parable to those of the individual single-2k models, but 
with a total mass M :^ 5 x lO^M© and a half-mass radius 
Th — 0.3 pc. The resulting half-mass relaxation time, trh, 
is ~ 4.4 Myr, which enable s these clusters to experience 
collapse within ^ 1 Myr (jPortegies Zwart fc McMillanI 
[2(K)2l : [Giirkan et al.l[2Q0l . 

The models, single-2k and single- w6, have a higher ini- 
tial central density than those of observed star clusters. 
This choice of initial density is motivated by our require- 
ment that the cluster should experience core collapse 
within ~ 1 Myr, which is required for the clusters to 
effectively to produce runaway stars. However, the core 
density drops and becomes comparable to the observed 
ones after the core collapse due to the scattering of stars 
in t he cluster center. We discuss more details in section 

We additionally performed simulations of a single 
model with the same parameters as single-w6 but with 
a smaller lower-mass limit of the IMF. For our standard 
models, we adopted 1 Mq as a lower-mass limit, which is 
more massive than observed ones. We adopted O.47M0 
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as a lower-mass limit for this model (single-lm), and as a 
consequence the mean mass of single-lm is half of single- 
w6. 

Furthermore, we performed two additional series of cal- 
culations of single clusters with the same mass and IMF 
as for the single-w6 model. One of these models has 
Wo = 8 and half-mass radius — 1.4 pc; we call this 
model single- w8. This model initially has a density pro- 
file similar to R136. The core density for this model is 
1.6 X 10^ M0pc~^ and its relaxation time is ~ 43 Myr; 
this cluster is therefore not expected to experience core 
collapse within ^ .p^^f (jPortegies Zwart fc McMillanI 
I2QQ2I : iGiirkan et al.l l2QQ4f). The other model is with pri- 
mordial binaries; which we name single-pb. We assigned 
the most massive 3% of stars ( ^ SM©) as equal-mass 
primordial binaries with the binding energy of 3 kT, 
which corresponds to an orbital period of ~ 3 x 10^- 
10^ days. Here kT is the fundamental unit of kinetic 
energy and 1.5 kT corresponds to the average kinetic 
energy of stars in the cluster before we add binaries. Bi- 
naries with this binding energy have intermediate hard- 
ness and therefore interact efficiently with single stars 
and with each other (Tanikawa & Fukushige 2009) and 
therefore they are expected to produce runaway stars, 
while binaries with shorter orbital periods are dynami- 
cally inactive due to their small cross sections (FPZ2011). 
We choose equal-mass binaries because massive stars ob- 
servatio nally prefer massive co mpanions (Sana & Evans 
[20111 : IZinnecker Yorkd 120071 ) . Ah binaries are initial- 
ized with circular orbits and random orbital phase and 
inclination. In table [2] we summarize the initial condi- 
tions for these models. 

During our simulations, we applied mass loss due to 
stellar winds only for stars more than lOOM©. Initially 
our simulations do not contain such massive stars, but 
they form occasionally due to collisions, for which we 
adopted the sticky-sphere approach. For the rate of mass 
loss by the stellar wind, for stars with m > 100M(7), we 
adopted m = 5.7 x lO~^(m/M0) (Moyr"^) (|Fuiii et al.l 
[2009). The stellar radii were adopted to have zero- age 
main-sequence radii for solar metalicity (jHurlev et al.l 
[2000). 

Each simulation was performed up to an age of 3 Myr, 
which is the time for the ffist supernova explosion. The 
equations of motion were integrated using a sixth-order 
Hermite scheme with individual timesteps with an accu- 
racy parameter 77 = 0.15 - 0.3 (Nitadori & Makino 2008). 
We adopted the accuracy parameter to balance speed and 
accuracy. Our code does not include any special treat- 
ment for hard binaries. The sixth-order Hermite scheme, 
however, can deal with hard binaries which form in our 
simulations. We tested our code integrating hard bina- 
ries with an eccentricity of > 0.9. Figure [T] shows an 
orbital separation of a binary. The masses of the binary 
components are 100 and 150 M©, respectively, and their 
orbital period, semi-major axis, and eccentricity are 72 
year, 109 AU, and 0.96. This binary has the same pa- 
rameters as one of hard binaries which actually formed 
in our simulations. The energy and angular momentum 
errors maintain < 10~^ and < 2 x 10~^ in 0.1 Myr. Even 
if the energy error increases linearly, the relative error in 
the energy is conserved to less than 0.1% over the 3 Myr 
of the simulation: our code treats close stellar encoun- 
ters down to ^ 1 AU accurately. Stars that approach 



each other within 1 AU generally participate in a colli- 
sion and therefore such close encounters do not need to 
be resolved. For each simulation we conffimed that the 
energy error of clusters is below 0.1% of the initial total 
energy of the clusters. 

3. COMPARISON BETWEEN SIMULATIONS AND 
OBSERVATIONS 

After 3 Myr, the sub-clusters merged to one cluster for 
all merger simulations. We present the time evolution 
of model m2k8r6 until 3 Myr in Figure [21 The sub- 
clusters hierarchically merge and finally form one star 
cluster with a number of escaping stars. The merger 
remnants are almost spherical and look very similar to 
initially spherical single clusters. However, the distribu- 
tion of stars in the merged clusters are different from 
that of clusters that were formed as one. In this section 
we demonstrate that the merger scenario is preferable for 
the formation of young dense clusters. 

3.1. Cumulative number of massive stars 

The projected cumulative distributions of massive 
stars as a function of the distance from the cluster center 
at a cluster age of 3 Myr are presented in Figu re [3l We 
adopted the density center (jCasertano fc Hut[ [l985) as 
the cluster center. The black dashed and dotted curves 
show the results for the single clusters, whereas the col- 
ored curves show the results for the merger simulations. 
We plot massive stars with > 3OM0 for the simulations. 
For comparison, we overlay the observed distribution for 
massive stars Wolf-Rayet (WR), spectral type O, and 
some e arly B stars from R136 (le ft) and from NGC 3603 
(right) ('Crowther fc Dessart' 1998, thick black curve with 
crosses). Although the mass range we adopted for the 
simulations is higher than that of the observation, the 
fraction of plotted stars to the total mass of the cluster 
is similar for both the simulation and the observation. If 
we assume that a Salpeter mass function between 0.3-100 
M0, the fraction of 0-stars (> I6M0) is 13.5%, which 
corresponds to the fraction of stars with > 3OM0 assum- 
ing a Salpeter mass function between 1-100 Mq. 

In the left panel of Figure [3l we find that the merger 
models and the single cluster with a relatively low cen- 
tral density (single- w8) are consistent with R136. We 
quantify the spatial distribution of the massive stars us- 
ing the radius of the inner tenth star, rio- These values 
are 0.11 ± 0.02, 0.087, and 0.15 pc for merger models, 
single- w8, and R136, respectively. Model single- w8 does 
not experience core collapse sufficiently quickly to ex- 
plain the high degree of mass segregation inside the star 
cluster, nor does this model produce the observed num- 
ber of runaway stars (the detail is in section 3.2). On 
the other hand, single models with a high central den- 
sity (single-w6) and with primordial binaries (single-pb) 
can produce enough runaway stars, but produce too a 
concentrated distributions of massive stars. Their values 
of rio are 0.034 ± 0.005 and 0.025 pc for single-w6 and 
single-pb, respectively, and these values are much smaller 
than that of R136. We also find the agreement between 
merger models and observation for NGC 3603 (see right 
panel of Figure [3]). 

We also see in the left panel of Figure [3] that a wider 
initial distribution of sub-clusters results in a relatively 
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TABLE 1 
Single cluster models. 



Model 


N 


M (Mo) 


Wo 


rh (pc) 


Pc (M0PC-3) 


cr (km/s) 


trh (Myr) 


rrimin (Mq) 


A^run 


single-2k 


2048 


6.3 X 103 


2 


0.097 


2.2 X 10^ 


11 


0.37 


1.0 


7 


single-w6 


16384 


5.1 X 10^ 


6 


0.32 


1.7 X 10^ 


17 


4.4 


1.0 


6 


single-lm 


32768 


5.1 X 10^ 


6 


0.32 


1.7 X 10^ 


17 


8.3 


0.47 


3 


single-w8 


16384 


5.1 X 10^ 


8 


1.4 


1.6 X 10^ 


6 


43 


1.0 


1 


single-pb 


16875 


5.1 X 10^ 


6 


0.32 


1.7 X 10^ 


17 


4.4 


1.0 


1 



TABLE 2 
Merger models. 



Model A^ei^ 


rmax (pc)"^ 


ATrun" 


m2k8rl 8 


1 


1 


m2k8r3 8 


3 


1 


m2k8r5 8 


5 


2 


m2k8r6 8 


6 


2 


m2k4r3 4 


3 


3 



Note. — The models are named according to the following rules; "m" indicates merger models, the name of models for the sub-clusters, 
the number of sub-clusters, and the value of rmax- For example, m2k8rl stands for a merger model of eight single-2k clusters which are 
located within rmax = 1 pc. 



^The number of sub-clusters. 

^The radius in which sub-clusters initially distribute. 
*^The number of runs which we performed. 




t (year) t (year) 

Fig. 1. — Orbital separation of a hard binary as a function of time (right for 0.1 Myr and left: for one orbital period). The masses of 
the binary components are 100 and 150 Mq, respectively. The orbital period, semi-major axis, and eccentricity are 72 year, 109 AU, and 
0.966, respectively. In right panel, crosses indicate the time steps adopted by our timestep criterion ([Nitadori &: Makindi2008i ') with an 
accuracy parameter of 0.15. The minimum and maximum stepsizes in this period are 0.73 day and 1.0 year, respectively. 



wider distribution of massive stars after the mergers, i.e., 
the value of rio is larger with the larger rmax- While 
rio = 0.096 ± 0.002 pc for the models with rmax = 1-3 
pc, rio = 0.12 ± 0.03 pc for the models with rmax = 5-6 
pc. The wider models can reproduce the distribution of 
massive stars in R136 better. 

We find that the degree of dynamical evolution in each 
of the sub-clusters before they merge makes this differ- 
ence. The sub-clusters merge after they experience core 
collapse for the majority of our simulations except model 
m2k8rl. In the latter case, the time scales for merger and 
core collapse are comparable. One would naively expect 
that the more dynamically developed sub-clusters will 
cause the merger product to be more centrally concen- 
trated also, but that is only partially true. The dynami- 



cally evolved sub-clusters results in the lower concentra- 
tion because of their shorter relaxation time. Figure |4] 
shows the time evolution of the core density for single 
clusters. The isolated sub-clusters (single-2k) collapse 
at ~ 0.3 Myr and then the core density decreases dra- 
matically due to the scattering of stars by hard binaries 
formed in a sub-clusters before they merged. This lower 
core density after the core collapse induces the lower con- 
centration of the merger remnants. Contrary to this, ini- 
tially massive single clusters maintain their high densities 
much longer. 

The wider initial distribution of sub-clusters provides 
sufficient time for them to experience core collapse and 
produce a BB, which ejects massive stars from each sub- 
cluster (FPZ2011). The post-collapse sub-clusters have 
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a lower density than those in a state of core collapse. Af- 
ter the sub-clusters merge, however, most of the BBs are 
broken or scattered up d ue to binary-binary interactions 
(jGualandris et al.ll2QQ4f ). and the production of massive 
runaways is finally driven by one single BB. The remain- 
ing cluster starts to evolve as an ordinary single cluster. 
This explains the wider radial profile of the more massive 
stars. We schematically represent this process in Figure 

El 

In Figure [6] we present the density profile of the sim- 
ulated clusters at 3 Myr. The central density of the 
merger products has already been depleted except model 
m2k8rl. These merger models and the single model 
single-w8 show similar distributions of massive stars and 
density profiles. Only the most concentrated merger 
model, m2k8rl, has parameters similar to a single dense 
model, single- w6. The central density (r ^ 0.2 pc) 
of the merger cases with rmax =5-6 pc is similar to 
the observed central dens ity of R136, 3-5x10^ MqPc"^ 
(jMackey fc GilmQr3l2QQ3f ) , and the central density of the 
four-merger cases is comparable to th e core density of 
NGC 3603, > 6 X 10"^ Mopc-^ (Harayama et al."2008'). 

From the too concentrated distribution of massive 
stars, we can exclude the single high density models 
(single- w6 and single-pb) as viable models of R136. If 
we see only the radial distribution in the cluster, how- 
ever, the difference between the merger and single cases 
is not clear. 

3.2. Mass distribution of runaway stars 

Another way to evaluate the degree of the dynamical 
evolution of young star clusters is finding runaway stars 
around the cluster. We demonstrated in FPZ2011 that 
the fraction of runaway stars around R136 and also their 
mass distribution are consistent with those scattered by 
a binary in the core of the host cluster. 

In FPZ2011, we discussed how a single cluster can ex- 
plain the fraction of runaway stars around R136. The 
degree of mass segregation in our model single-w6 is too 
strong compared with the observed distribution of mas- 
sive stars (see Figure [3]). On the other hand, if we as- 
sume a single cluster with a distribution of massive stars 
(single-w8) similar to R136, no runaway stars are pro- 
duced because the relaxation time of this model {trh = 44 
Myr) is too long to experience core collapse within 3 Myr 
(see Figure 3]). 

We conclude that both the observed distributions of 
massive stars and runaway stars of R136 can be repro- 
duced if the star cluster formed through the hierarchi- 
cal merger of ~ 8 small sub-clusters. In Figure [3 we 
show that /run, the fraction of runaway stars that es- 
cape their host cluster with > 30km/s, agrees well with 
the WR and spectral type O- and early B-stars observed 
around R136 (square). The observed /run of early B- 
type stars (8-16 Mq) is a few times smaller than that 
expected from our si mulations. Assuming that th e total 
mass of 6.0 x lO'^M© (jPortegies Zwart et al.ll2010f ) with a 
Salpeter mass function between 1 and 100 Mq and that 
the runaway fraction obtained from our 8-merger mod- 
els, we predict that there are 7±2 early-B runaways still 
in the vicinity of R136. We also expect that ~ 6 late-0 
(16-32 Mq) and - 11 early-O/WR (> 32M0) runaway 
stars exist around R136. 

The lower-mass limit of the IMF does not affect on the 



fraction of runaway stars (see the right panel of Figure [7]) 
because the three-body scattering by a binary in the clus- 
ter core is driven by massive stars gathering in the cluster 
core due to the mass segregation. Even if we change the 
lower- mass limit of the IMF, low-mass stars do not in- 
teract much frequently with the binary in the core. In 
Figure [8] we show the kinetic energy of runaway stars 
as a cumulative function of mass, which corresponds to 
the energy extracted from the binary by three-body scat- 
terings. The contribution from low-mass stars is much 
smaller than that from massive stars. The core-collapse 
time is not sensitive to the lower- mass limit either, al- 
though the number of particles increases with a smaller 
lower-mass limit. From the results of A/'-body simula- 
tions, the core-collapse time, tec, follows an equation such 
as tec {mma^/ {m))~^-^trh (Giirkan et al. 2004), where 
mmax and (m) are lower- mass limit and mean mass, re- 
spectively. Since (m) ex 1/A^ and the relaxation time 
is roughly proportional to A/", the core-collapse time de- 
pends on A""^-^. The lower-mass limit of the IMF, how- 
ever, changes the mass where the fraction of runaway 
stars starts increase because only stars which are more 
mas sive tha n t his l imit mass can sink to the cluster cen- 
ter (jFujii et al.ll2008f ). A smaller value of the lo wer-mass 
limit decrease this mass (see Figure [71 and Ba nerjee et al.l 
I2OI2L in which they assumed an IMF between 0.3 and 150 
Mq and obtained very similar results to ours.). The limit 
mass is 13 and 9 Mq for our models with a lower-mass 
limit of 1 and 0.47 M©, respectively. 

The consistency between the observed number of run- 
aways and the density profile for R136 gives us confidence 
in our prediction of a relative fraction of runaway stars 
around NGC 3603 and other clusters. As we discussed 
before, the fraction of runaway stars does not depends on 
the lower-mass limit of the IMF. Therefore, we adopted 
both 0.3 and 1 Mq as the lower-mass limit and assumed 
the fraction of runaway stars obtained from our simula- 
tions a the lower-mass limit of 1 Mq. Since the number 
of massive stars assuming a lower-mass limit of 0.3 Mq 
is ~60% of that in the case of 1 M©, the expected num- 
ber of runaway stars also decrease ~60%. Adopting a 
Salpeter mass function between 1 Mcd and 100 M^T) and a 
total cluster mass of - l.TxlO'^M© (|R.ochau et al.ll2010f ). 
we predict that — 7 0/WR stars (> I6M0) have been 
ejected from NGC 3603. If we assume a lower- mass limit 
of 0.3 M0, we predict — 4 0/WR runaway stars. We 
summarize the expected number of OB runaway stars 
around young star clusters in table [3l We assumed 8- 
merger models for R136 and Westerlund 1 and 4- merger 
models for NGC 3603 and Westerlund 2. 

3.3. Single cluster with primordial binaries 

Since the binary fraction of massive stars in star clus- 
ters is very high (20-80%) (|Zinnecker fc Yorkd 120071 . 
and references therein), and binary-binary interac- 
tions are common and can e ffectively produce run- 
aways (Leonard & Duncan 1988: iClarke fc PringldfTool : 
Gualandris et al. 2004), a fraction of primordial binaries 
and their distribution in binding energies are expected 
to affect the fraction of runaway stars. We performed 
a simulation for a single cluster model with primordial 
binaries (single-pb) in order to quantify the effect of bi- 
naries. The radial distribution of massive stars (Figure 
[3]) and the density profile (Figure [6]) are very similar to 
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the single case with the same initial density distribution 
(single- w6). The primordial binaries do not solve the 
problems that clusters require too high a core density 
and too concentrated a distribution of massive stars to 
form a sufficient amount of runaway stars; primordial bi- 
naries therefore do not help in reproducing the spatial 
distribution of the stars in R136. 

The density profile of the cluster is not affected by the 
presence of primordial binaries, but the fraction of the 
runaway stars is higher compared to those produced in 



single clusters, and similar to those produced in merged 
cluster (see Figure [7j). This fraction is cons i stent with 
the simulations performed by iBanerjee et al.l (|2Q12l ) . 

3.4. Comparison with Westerlund 1 

Westerlund 1 is another young massive star cluster in 
the Milky Way with a mass comparable to R136, but 
with an age of 3-6 Myr ()Clark et al.ll2QQ5l : iGennaro et al.l 
I2Q11I ). The relaxation time estimated from its current 
mass and size is ^ 130 Myr, b ut the cluster is mass- 
segregated (Gennaro et al.|[2QlT| ). Furthermore, Wester- 
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Fig. 3. — Projected cumulative number distribution of the massive (> SOMq for simulations) stars for R136 (left) and NGC 3603 (right) 
models. Col ored curves sh ow merger models and black dashed, dash-dotted, and dotted curves show single models. Black crosses show the 
observation (jCrowther D essart 1998). 



TABLE 3 
Expected runaway stars. 



Name 


Mel (Mo) 


Ref. 


Model 


mmin(M0) 


early-B 


late-O 


early-O/WR 


R136 


6.0 X 10^ 


1 


8-merger 


1 


6.7 ±2.4 


6.4 ±3.0 


11 ±6 










0.3 


4.9 ±1.5 


3.9 ±1.8 


6.5 ±3.5 


NGC 3603 


1.7 X 10^ 


2 


4-merger 


1 


2.9 ±2.2 


2.1 ±0.6 


4.8 ±1.5 










0.3 


1.8 ±1.3 


1.3 ±0.3 


2.9 ±0.9 


Westerlund 1 


6.3 X 10^ 


1 


8-merger 


1 


7.0 ±2.5 


6.8 ±3.2 


11 ±6 










0.3 


4.2 ±1.5 


4.1 ±1.8 


6.8 ±3.7 


Westerlund 2 


10^ 


1 


4-merger 


1 


1.7 ± 1.3 


1.2 ±0.3 


2.8 ±0.9 










0.3 


1.0 ±0.8 


0.8 ±0.2 


1.7 ±0.8 



Note. — References: 1. IPortegies Zwart et al.) (|20Tol ). 2. IRochau et"all (|20Tol V The lower -mass limit is indicated as rrimy 




t (Myr) 

Fig. 4. — Time evolution of core density. These are averaged on 
all runs for models single-w6 and single-2k. 



lund 1 is highly asymmetric, which could indicate that 
it had experienced mergers. Unfortunately, there is no 
data of runaway stars around this cluster, and little data 
is available about its internal kinematics. However, the 
effective cumulative radius: 



reffM = \ , (1) 



has been observed (jGennaro et al.|[2QTT| ). Here, ri is the 
distance from the cluster center in projection for i-th 
massive stars with rui. This is the geometric- averaged 
distance from the cluster center for stars with mass 
rrii > m. We calculated reff(m) for our simulations in 
order to compare it with the observation. We took into 
account stars in a radius of 3pc to remove the effect of 
runaway stars. This radius is the same as the observed 
region, ~ 2.3 pc. The results are shown in Figure [9l The 
smaller effective radius for more massive stars indicates 
mass segregation. With simulation m2k8r3, which has a 
density profile which is consistent to the observed cluster 
R136, we quantify the observation to reff ~ 0.8 pc at the 
minimum radius at massive end (~ 2OM0) and rgff ~ 1.5 
pc for all stars dow n to 2.5 (see also Figure 11 of 
iGennaro et al.|[2QTTh . The effective radius for stars with 
^ 3OM0 is larger than than that at around 2O-3OM0 in 
both the simulation and the observation. In the simula- 
tion the larger radius comes from some off-center mas- 
sive stars, which are scattered from the cluster center. In 
contrast to the merger models, the single model single- 
w8, which initially has properties similar to of that of 
the current Westerlund 1, has not have sufficient time 
to evolve dynamically. From these results, we consider 
that Westerlund 1 is also a young massive cluster which 
is forming via hierarchical merger. Since Westerlund 1 
(3-6 Myr) is slightly older than R136 (- 3Myr) but is 
still undergoing mergers, we expect that the sub-clusters 
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Fig. 5. — Formation of a star cluster through mergers. Sub-clusters form in a region of several pc (left top). They dynamically evolve 
and experience core collapse within 1 Myr. A BB form in each cluster and produce runaway stars by three-body scattering (right top). 
The sub-clusters merge hierarchically and the hardest BBs break or kick out the other BBs (bottom). 




of Westerlund 1 were born in a larger region than those 
of R136. We expect that based on the simulation, the 
numbers of runaway stars of Westerlund 1 are 7.0 ± 2.5 
(4.2 ± 1.5), 6.8 ± 3.2 (4.1 ± 1.8), and 11 ± 6 (6.8 ± 3.7) 
for early-B, late-0, and early-O/WR stars assuming a 
Salpeter IMF with 1-100 (O.3-IOOM0), respectively 
(see also table [3j). 



4. DISTRIBUTION OF COLLIDED STARS 



Recent observations have identified several very mas- 
sive stars with an estimated zero-age mai n-sequence mass 
of ^ I5OM0 in R136 and NGC 3603 (Crowth er etHI 
|2010[ ). Because collisions between stars easily happen in 
a dense cluster (jPortegies Zwart et al.l ll999). such very 
massive stars are a natural by-product of the dynamical 
evolution of a dense star cluster. 

In Fig. [To] we present the radial distribution and the 
mass distribution of stars that experienced collisions in 
our simulations. We adopted for adopting the maximum 
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m(M0) m{MQ) 

Fig. 7. — Fraction of runaway stars (> 30 km/s) to whole cluster as a function of mass at 3 Myr. Color curves show the results 
of our simulations. Left and right panels show comparisons between single and merger models and among single models, respecti vely. 
Square and triangle symbols show the runaway fraction of R136 and the field (from FPZ2011; Gvaramadze et al. 2010; Bestenlehner e t al.l 
[201 li: t Andersen et al.ij2009i: iZinnecker &: Yorkei »2007|) . S t ar shows a fraction of is olated WR and early O-stars around Westerlund 2 
(ITsuiimoto et al.ll2007l : IRauw et al.ll2007l : INaze et al.ll2008l : IRoman-Lopes et al.|[20Tir) . 



^ 400 




Fig. 8. — Kinetic energy of runaway stars as a cumulative func- 
tion of mass. The energy is normalized by kT of each cluster. Note 
that the value of kT for model single-w6 is twice as large as that 
of single-lm because of the difference of their mean masses. 



2.0 



1.5 



single-w6 
single-w8 
single-pb 



m2k8r3 

■ - m2k8r5-1 
+ Westerlund 2 




10 

m(M0) 

Fig. 9. — Effective cumulative radius for stars in 3 pc (in 
projection) at 3.5 M yr. Data points of Westerlund 2 are from 
IGennaro etaLl (pOlli ). 



masses a star reaches during its hfetime rather than the 
mass at the end of the simulation, because the latter 
strongly depends on the mass-loss rate. For comparison, 
we also indicate the very massive stars observed in R136 
and NGC 3603 (iCrowther et alJl2QlQl : lBestenlehner et all 
Hoil). In the left panel of Figure [lOl we can see the effect 
of the mass segregation within ~ 1 pc, while we find stars 
escaping at distances greater than 10 pc. The escaping 
massive stars formed in the merger models with rmax = 
5-6 pc. The maximum mass that was produced in these 
simulations was around 24OM0, and three of these stars 
are ejected from the cluster as single runaways. Their 
velocities are similar, ~ 60 km/s in three dimensions, 
which corresponds to ~ 35 km/s in one dimension. This 
mass and velocity are very similar to those of VFTS 682, 
which is located at ^ 30 pc from R136 in projection and 
seems to be escaping f rom R136 with one-dime nsional 
velocity of - 30 km/s CBestenl ehner et al.l[2QTTh . With 
a current mass of VFTS 682 is ^ 15OM0 its initial mass 
would < 2IOM0 ([Bestenlehner et al.l[20TTh . 

Such very massive stars have been found also inside of 
R136. There are four massive star s with in itial masses 
of 165-320 Mq in R136 (|Crowther eHIIImO j. Three of 
four massive stars in R136 are located in the inner most 
0.1 pc and the other one is ~ Ipc (in projection) from 
the cluster center. In our simulations, we found 3-4 very 
massive stars in the cluster center in the case of merged 
clusters with rmax = 5-6 pc. Model m2k8r6-l evolves 
towards a distribution very similar to R136. In Figure 
[TH we present the radial distribution of merger products 
with > IOOMq. 

We compared the distribution of massive stars in NGC 
3603 with the simulations. There are four very mas- 
sive stars with initial masses of 105-170 Mq in NGC 
3603, and two of them (Ala and Alb) form a binary 
(ICrowther et all 120101 ). We found 2-3 collided stars 
which exceed 100 Mq in each run. Some of the stellar 
merger products form binaries. We find a massive bi- 
nary composed of stars that experienced a collision in 
model m2k4r3-3 and one composed of a collided star 
and a ~ lOOM© normal star in model m2k4r3-l. We 
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speculate, based on our simulation results, that massive 
(> IOOM0) star could have been ejected from NGC 3603 
(see right panels of figures [10] and [11]), and could still 
roam the LMC. 

The merger simulations that started with a relatively 
large rmax have mass and density profiles that are con- 
sistent with the observations, but the single clusters and 
merger simulations with small rmax are unable to repro- 
duce the observations. In the latter cases, one very mas- 
sive star ^ 4OOM0 forms in the cluster as a result of 
runaway collisions (see Figure [TO]). However, there is no 
observational evidence of such a massive star in R136 or 
NGC 3603. 

In the simulations with primordial binaries (single- 
pb), no massive runaway star which experienced colli- 
sions formed, whereas Ban eriee et al.l ^2012) found some 
collided runaway stars in their A^-body simulation mod- 
els of R136 for a single cluster with primordial bina- 
ries. In our simulation, most of the massive stars col- 
lided into the most massive star via runaway collisions 
as is the c ase without primordial b inarie s, which is consis- 
tent with iPortegies Zwart et al.l (|1999[ ) . We argue that 
th e difference bet ween our results and those presented 
in lBanerjee et al.] (20T2) originate from the initial condi- 
tions. The major differences are the orbital period distri- 
bution of primordial binaries and their assumption that 
all binaries are born in the cluster core. We adopted the 
same binding energy for all binaries (orbital periods for 
massive binaries are ~ 3 x 10^-10^ days) and distribute 
them in the same w ay in the cluster as th e single stars. 
On the other hand, iBanerjee et al.l (|2012[ ) chose shorter 
periods of 0.5 < log^^g (^/^^Y^) < 4. They employed ini- 
tial mass segregation; they initialized the heaviest stars 
in the cluster core. The combination of relatively tight 
binaries and a high degree of mass segregation stimu- 
late strong dynamical encounters which may in their case 
have resulted in the stimulated formation of a massive 
runaway in the simulations with primordial binaries. 

In Figure [11] we show that about half the stars that 
experienced one or more collisions are accompanied by 
another star. In simulations with merging clusters, one 
BB forms in each sub-cluster before it merges. The BBs 
in the sub-cluster centers participate in colliding encoun- 
ters due to which they become more massive. After the 
sub-clusters merge, one of the BBs can be dynamically 
scattered and escapes from the cluster. We discuss bina- 
ries in section ??. 

Stars that e xperience a collision are candidates for 
rapid rotation (jFrver fc He gerl [20051 ). Several rapidly ro- 
tating massive stars have been found in recent observa- 
tions, and possibly this rotation is the result of the earlier 
collision history of the star. The star VFTS 682, has a 
rotation velocity > 200 km/s (jBestenlehner et al.l 1201 if ), 
and the > 2OM0 VFTS 102, which is a runaway with 
~ 40 km/s located in the 30 Doradus refflon, is also a 
rapid rotator (> 500 km/s) (jPufton et al.l 1201 ID . The 
rapid rotation of these stars could be explained by an 
earlier episode of mass transfer in a binary system or by 
collision with another star. Either of these probably oc- 
curred in the cluster center, from which it can easily have 
been ejected. In our simulations we find several rotating 
stars inside the model that reproduces the observations 
of R136 and a few are ejected as runaways (see figures 



[IO]and[IID. 

5. DISTRIBUTION OF BINARIES 

In Figure [12] we show the radial distribution of hard 
binaries (the binding energy is more than IkT) at 3 Myr. 
In the simulations with primordial binaries, we find 306 
such binaries at an age of 3 Myr, which is ~ 60% of 
the primordial binary population. The majority of these 
binaries are primordial. 

In the case of a single cluster with/without primordial 
binaries and merged clusters with rmax = 1-3 pc, there 
is generally one BB in the cluster center. In the case of 
a merged cluster with rmax = 5-6 pc, however, we found 
3-7 BBs. Furthermore, we found two escaping binaries 
in model m2k8r6-l. They are located at 38 and 27 pc 
from the cluster center and their velocities are 24 km/s 
and 16 km/s, which are slightly slower than the velocity 
to be defined as runaway stars, but these are unbound. 
We also found off-center binaries in two of our merger 
models (models m2k8r5-l and m2k8r6-2). They are lo- 
cated at 36 and 8 pc from the cluster center, but their 
velocity is ~ 10 km/s and they are bound to the clus- 
ter. These binaries are remnants of the BBs that formed 
in sub-clusters and were kicked out after their host clus- 
ters merged. Therefore, they are very massive with a 
maximum mass ~ 24OM0. Figure [13] shows the orbital 
periods and the total current masses of the binaries ob- 
tained from the simulations. The orbital period of the 
massive binaries located farther than 10 pc (marked by 
circle in the figure) is ~ 500 days. R145, which is a mas- 
sive binary (300 + 125 Mq) located at ~ 20 pc in pro- 
ject ion_from_^136and_of which orbital period is 158. 8 
days ( Schnurr et al.l 120091 : iFitzpatrick fc Savagd [19841 ). 
has similar characteristics to these escaping binaries. 

6. SUMMARY 

We performed simulations of merging and single star 
clusters and compared the radial distribution of massive 
stars and the mass distribution of runaway stars to obser- 
vations. We found that clusters that assembled from the 
hierarchical merging of multiple sub-clusters can explain 
the characteristics (mass distribution, mass segregation, 
runaway stars, massive stars, and massive binaries) of 
young dense clusters such as R136 in the LMC and NGC 
3603, Wester lund 1, and 2 in the Milky Way. 

We have set-up our simulations in such a way that the 
sub-clusters have a sufficiently short relaxation time that 
they experience core collapse before they merge and as a 
consequence the large merged cluster is in a state of core 
collapse when it forms. This results in the mass segre- 
gation of the merger remnants. Each sub-cluster experi- 
ences core collapse which leads to a collision runaway and 
the consequent production of a very massive (> IOOM0) 
star in the process. After the clusters have merged the 
BBs gather in the center of the merged cluster. The ma- 
jority of them are ionized or ejected from the cluster out 
by binary-binary, in fact BB-BB, interactions. These 
processes lead to the formation of very massive stars, 
consistent with those observed in R136 and NGC 3603, 
massive runaway stars such as VFTS 682, and massive 
binaries in the periphery of young cluster like R145. The 
post-core-collapse evolution of the sub-clusters drives the 
expansion of the core by the ejection of a large number 
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Fig. 10. — Radial and mass distribution of collided stars at 3 Myr for 8- merger models (left) and 2.5 Myr for 4-merger models (right). 
Circles show merged cluster cases, whereas squares show a single case with primordial binaries. The cross shows the averaged value of single 
cluster models without primordial binaries (single-w6). In the case of single cluster s, only one runaway co llision star g rows in a cluster. 
Black triangles show very massive stars (> lOOM©) found in R136 and NGC 3603 f Crowther et al.l 12010 : IBestenlehner et al.ii2011l ). The 
right top green point in the right panel is a binary which have very similar masses. 
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Fig. 11. — Cumulative radial distribution of stars which experienced collisions and with > lOOM© at 3 Myr for R136 models (left) and 
2.5 Myr for NGC 3603 models (right). Crosses indicate binaries whose companions are also collided stars (two crosses at the same radius 
are a pair). Surrounded circles indicate binaries whose companions are not collided stars. 



of runaway stars. The distribution of massive stars and 
the core density of the star clusters which formed via hi- 
erarchical mergers are consistent with those of R136 and 
NGC 3603. 

The initial conditions of our simulations are con- 
structed such that core collapse occurs at an age of 
~0.3Myr. This fine tuning is required in order to ex- 
plain the coexistence of a rich population of runaways, 
the presence of a BB in the cluster center (or outside) 
and the relatively low density of the observed clusters. 
Observing a star cluster in a state of cluster collapse at 
an age ^ 1 Myr would provide a confirmation for forma- 
tion model for massive clusters and the co-production of 
runaway stars and extremely massive single stars. 

If we adopt single clusters with initial conditions for 
which the distribution of the massive stars is consis- 
tent with observations, our simulated clusters are insuffi- 
ciently dense to experience core collapse within a few Myr 
and to form runaway stars within the age of the observed 



clusters. On the other hand, initial single clusters that 
are sufficiently dense to experience core collapse within 
3 Myr can produce a sufficiently large number of run- 
aways compared to the observed number found around 
R136, but they are considerably more concentrated than 
the observed clusters. Single clusters with primordial bi- 
naries mediates the number of runaway stars, but also 
lead to a too concentrated density proffle at later time. 

Westerlund 1 and 2 are also candidates of merged 
clusters, in particular because they appear to be 
rather clumpy and t hey are strongly m ass segregated 
(jGennaro et 8011201 II : lAscenso et al.l 1200 7). The degree 
of mass segregation in our merger simulations is con- 
sistent with the observed cumulative effective radius of 
Westerlund 1 . The fraction of runaway stars in the vicin- 
ity of Westerlund 2 is consistent with the result of our 
simulations in which multiple clusters merge to one. 

Our results support the claim that some young dense 
clusters, such as R136, NGC 3603, Westerlund 1, and 2, 
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r (pc) 

Fig. 12. — Cumulative radial distribution of hard binaries at 3 
Myr. Colored points show merger models, and the black curve 
shows a single model with primordial binaries. 



formed via mergers of smaller clusters. We further argue 
that the small sub-clusters which are best suited to re- 
produce the observed characteristics of the young dense 
clusters are also the building blocks for star clusters else- 
where in the Milky Way, and are therefore responsible for 
the production of O and B runaway stars in the Galactic 
disk. The typical scale of the cluster is M ~ GOOOM©, 
and the fraction of runaway stars among 0-type stars 
is about 10-20%, and drops to a few percent for B-type 
stars (FPZ2011). Clusters more massive than those stud- 
ied here may also have formed through the merging of 
multiple smaller sub-clusters. We argue that the funda- 
mental building-block cluster with which we are able to 
explain the density profile of observed star clusters and 
at the same time the rich population of massive early B 
and O stars in the vicinity has a mass of ^ lOOOOM©. 
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Fig. 13. — Orbital periods and total masses of binaries at 3 Myr. 
Points with circles indicate binaries which are located at > 5 pc 
from the cluster center. 
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